%=== Illustrate how changing the parameter c affects the relative benefit of entering each job type %
function [t_Uhat_N_compare,t_Uhat_O_compare,c_grid,num_c] = DO_c_effect(param)

    %----- Specify Parameters -----%
    nu = param(1);
    beta = param(2);            % Discount factor (factors in exit probability)

    alpha = param(3);           % Production function parameter
    d = param(4);               % Human capital depreciation rate
    ell = param(5);             % Matching function parameter
    lambda = param(6);          % Probability of search OTJ
    chi = param(7);             % Investment cost parameter
    b=param(8);                 % Unemployment benefit

    muA =param(9);              % Staffing agency fee
    %c =param(10);               % Cost of posting a vacancy
    delta =param(11);           % Separation rate
    h_sigma =param(12);         % Standard deviation of initial human capital distribution

    pr_phi(1) = param(13);      % Probability of entering the model with phi=1 (low fixed productivity type)
    pr_phi(2) = param(14);      % Probability of entering the model with phi=2 (high fixed productivity type)

    z(1) = param(15);           % Productivity parameter f(phi,h) = z(phi)h^(alpha)
    z(2) = param(16);

    muA_param=muA;
    gamma=2;
    vfi_toler = 0.0001;
    dist_toler = 0.0000001;
    
    % Specify grids %
    nphi = 2;                  % Number of possible phi values (fixed productivity types)
    hlb = 0.0;                 % Human capital grid lower bound 
    hub = 7.0;                 % Human capital grid upper bound 
    nh = 3000;                 % Number of possible human capital values
    h_grid = linspace(hlb,hub,nh); % Create human capital grid
    
    % ======================= Create discretized h distribution ==================== %
    h_mu = 1;                  % Mean of initial human capital distribution
    cdf_h = zeros(1,nh);       % CDF of initial h distribution
    pdf_h = zeros(1,nh);       % PDF of initial h distribution

    h_bin_size = (hub-hlb)/(nh-1);
    for ind_h = 1:nh
        if (ind_h ==1)
            bin_ub_now = h_grid(ind_h) + (h_bin_size)/2;
            cdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma);
            pdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma);
        else
            bin_ub_now = h_grid(ind_h) + (h_bin_size)/2;
            bin_ub_last = h_grid(ind_h-1) + (h_bin_size)/2;
            cdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma);
            pdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma) - normcdf(bin_ub_last,h_mu,h_sigma);
        end
    end


    num_c = 20;
    c_grid = linspace(0.1,1,num_c);
    t_Uhat_O_compare = zeros(nphi,nh,num_c);
    t_Uhat_N_compare = zeros(nphi,nh,num_c);
    for ind_c = 1:num_c
        c = c_grid(ind_c);

        %% ============== Value Function Iteration to Get Value and Policy Functions ================ %%
        %---- Prespecify Value Functions
        
        U = zeros(nphi,nh);                % Value of unemployment in production phase                     
        U_hat  = zeros(nphi,nh);           % Value of unemployment in search phase - when searching for optimal job type
        U_hat_update  = zeros(nphi,nh);
        
        VN = zeros(nphi,nh);               % Value of employment in non-outsourced (N) job in production phase
        VO = zeros(nphi,nh);               % Value of employment in outsourced (O) job in production phase
        VN_hat = zeros(nphi,nh);           % Value of employment in N job in search phase
        VO_hat = zeros(nphi,nh);           % Value of employment in O job in search phase
        VN_hat_update = zeros(nphi,nh);
        VO_hat_update = zeros(nphi,nh);
        
        A_staffing = zeros(nphi,nh);
        
        %---- Prespecify Value Functions
        opt_ind_hp_N = ones(nphi,nh);     % Optimal index for next period h value (h') if in non-outsourced job
        opt_ind_hp_O = ones(nphi,nh);     % Opitmal index for next period h value (h') if in outsourced job
        
        GU = zeros(nphi,nh);              % Indicates optimal job type to search for if unemployed - 1=non-outsourced (N) and 2=outsourced (O)
        GN = zeros(nphi,nh);              % Indicates optimal job type to search for if employed in N job - 1=non-outsourced and 2=outsourced
        GO = zeros(nphi,nh);              % Indicates optimal job type to search for if employed in O job - 1=non-outsourced and 2=outsourced
        
        opt_p_U = zeros(nphi,nh);         % Optimal job-finding rate associated with optimal search choice when unemployed
        opt_p_N = zeros(nphi,nh);         % Optimal job-finding rate associated with optimal search choice when employed in N job
        opt_p_O = zeros(nphi,nh);         % Optimal job-finding rate associated with optimal search choice when employed in O job
        
        opt_x_U = zeros(nphi,nh);         % x contract value associated with optimal search choice when unemployed
        opt_x_N = zeros(nphi,nh);         % x contract value associated with optimal search choice when employed in N job
        opt_x_O = zeros(nphi,nh);         % x contract value associated with optimal search choice when employed in O job
        
        opt_eta_U = zeros(nphi,nh);       % eta (surplus rate) associated with optimal search choice when unemployed
        opt_eta_N = zeros(nphi,nh);       % eta (surplus rate) associated with optimal search choice when employed in N job 
        opt_eta_O = zeros(nphi,nh);       % eta (surplus rate) associated with optimal search choice when employed in O job
        
        opt_theta_U = zeros(nphi,nh);     % theta (market tightness) associated with optimal search choice when unemployed
        opt_theta_N = zeros(nphi,nh);     % theta (market tightness) associated with optimal search choice when employed in N job 
        opt_theta_O = zeros(nphi,nh);     % theta (market tightness) associated with optimal search choice when employed in O job
        
        % Find index for next period human capital (h') if the agent only lets human capital depreciate
        ind_hp_min = ones(1,nh);          % Minimum index for h' (index for h' if human capital only depreciates)
        for ind_h = 1:nh
            h = h_grid(ind_h);
            error_min = 100;
            for ind_hp = 1:nh
                hp = h_grid(ind_hp);       % hp is the value of h'
                error = abs(hp - h*(1-d)); % Find index that minimizes (hp-h*(1-d))
                if (error<error_min)
                    error_min = error;
                    ind_hp_min(ind_h) = ind_hp;
                end
            end
        end
        
        %--------------- Begin Value Function Iteration ----------------%
        vfi_error = 1; % Set vfi_error > vfi_toler
        while (vfi_error > vfi_toler)
        
            for ind_phi = 1:nphi          % Loop over phi types
                for ind_h = 1:nh          % Loop over current h values
                    h = h_grid(ind_h);    % h = current human capital (h) value 
                    
                    b_effect = b*z(ind_phi)*(h^alpha);  % Unemployment value
                    c_effect = c;                                % Search cost
                    
                    %===== For section 5.4 =====% (all other policy changes only affect parameters stored in param)%
                    % ---------- Solve hc Investment Decision ----------%
                    max_VN = -9999;
                    max_VO = -9999;
        
                    for ind_hp = ind_hp_min(ind_h):nh      % Loop over possible next period hp (h') values
                        hp = h_grid(ind_hp);
        
                        temp_VN = z(ind_phi)*(h^alpha) - ((hp - h*(1-d))*chi)^gamma ...
                            + beta*(1-delta)*VN_hat(ind_phi,ind_hp) + beta*delta*U_hat(ind_phi,ind_hp);
        
                        temp_VO = (1-muA)*(z(ind_phi)*(h^alpha)) - ((hp - h*(1-d))*chi)^gamma ...
                            + beta*(1-delta)*VO_hat(ind_phi,ind_hp) + beta*delta*U_hat(ind_phi,ind_hp);
        
                        if (temp_VN>max_VN)                % Save ind_hp where N job employment value is highest
                            max_VN = temp_VN;
                            VN(ind_phi,ind_h) = temp_VN;
                            opt_ind_hp_N(ind_phi,ind_h) = ind_hp;
                        end
                        if (temp_VO>max_VO)                % Save ind_hp where O job employment value is highest
                            max_VO = temp_VO;
                            VO(ind_phi,ind_h) = temp_VO;
                            opt_ind_hp_O(ind_phi,ind_h) = ind_hp;
                        end
        
                    end % end loop over hp options
                    U(ind_phi,ind_h) = b_effect + beta*U_hat(ind_phi,ind_hp_min(ind_h));
                    
                    ind_hp = opt_ind_hp_O(ind_phi,ind_h);
                    A_staffing(ind_phi,ind_h) = muA*z(ind_phi)*(h^alpha) + (1-delta)*beta*(1-lambda*opt_p_O(ind_phi,ind_h))*A_staffing(ind_phi,ind_hp);
    
        
                    %---------------------------------------------------%
                    % ------- Solve Search Decision of Unemployed ----- %
                    xO = VO(ind_phi,ind_h) - c_effect;
                    etaO = (xO - U(ind_phi,ind_h))/(VO(ind_phi,ind_h) - U(ind_phi,ind_h));
                    t_Uhat_O = xO;
        
                    thetaN = (((VN(ind_phi,ind_h) - U(ind_phi,ind_h))/c_effect)^(ell/(1+ell)) - 1)^(1/ell);
                    if (~isreal(thetaN))
                        thetaN = 0;
                        pN = 0;
                        xN = 0;
                        etaN = 0;
                    else
                        pN = thetaN/((1+thetaN^ell)^(1/ell));
                        qN = pN/thetaN;
                        xN = VN(ind_phi,ind_h) - (c_effect/qN);
                        etaN = (xN - U(ind_phi,ind_h))/(VN(ind_phi,ind_h) - U(ind_phi,ind_h));
                    end
                    t_Uhat_N = pN*xN + (1-pN)*U(ind_phi,ind_h);
        
                    % Now Consider which job type's value is highest
                    if (t_Uhat_N>t_Uhat_O)
                        GU(ind_phi,ind_h) = 1;
                        opt_p_U(ind_phi,ind_h) = pN;
                        U_hat_update(ind_phi,ind_h) = t_Uhat_N;
                        opt_x_U(ind_phi,ind_h) = xN;
                        opt_eta_U(ind_phi,ind_h) = etaN;
                        opt_theta_U(ind_phi,ind_h) = thetaN;
                    elseif (t_Uhat_O>t_Uhat_N)
                        GU(ind_phi,ind_h) = 2;
                        opt_p_U(ind_phi,ind_h) = 1;
                        U_hat_update(ind_phi,ind_h) = t_Uhat_O;
                        opt_x_U(ind_phi,ind_h) = xO;
                        opt_eta_U(ind_phi,ind_h) = etaO;
                        opt_theta_U(ind_phi,ind_h) = 1.0;
                    else
                        GU(ind_phi,ind_h) = 3;
                        opt_p_U(ind_phi,ind_h) = pN;
                        U_hat_update(ind_phi,ind_h) = t_Uhat_N;
                        opt_x_U(ind_phi,ind_h) = xN;
                        opt_eta_U(ind_phi,ind_h) = etaN;
                        opt_theta_U(ind_phi,ind_h) = thetaN;
                    end
        
                    %---------------------------------------------------%
                    % ------- Solve Search Decision of Employed --------%
                    xO = VO(ind_phi,ind_h) - c_effect;
                    etaO = (xO - U(ind_phi,ind_h))/(VO(ind_phi,ind_h) - U(ind_phi,ind_h));
        
                    t_VNhat_O = lambda*xO + (1-lambda)*VN(ind_phi,ind_h);
                    t_VOhat_O = lambda*xO + (1-lambda)*VO(ind_phi,ind_h);
        
                    thetaN_N = (((VN(ind_phi,ind_h) - VN(ind_phi,ind_h))/c_effect)^(ell/(1+ell)) - 1)^(1/ell);
                    if (~isreal(thetaN_N))
                        thetaN_N = 0;
                        pN_N = 0;
                        xN_N = 0;
                        etaN_N = 0;
                    else
                        pN_N = thetaN_N/((1+thetaN_N^ell)^(1/ell));
                        qN_N = pN_N/thetaN_N;
                        xN_N = VN(ind_phi,ind_h) - (c_effect/qN_N);
                        etaN_N = (xN_N - U(ind_phi,ind_h))/(VN(ind_phi,ind_h) - U(ind_phi,ind_h));
                    end
                    t_VNhat_N = lambda*pN_N*xN_N + (1-lambda*pN_N)*VN(ind_phi,ind_h);
        
                    thetaO_N = (((VN(ind_phi,ind_h) - VO(ind_phi,ind_h))/c_effect)^(ell/(1+ell)) - 1)^(1/ell);
                    if (~isreal(thetaO_N))
                        thetaO_N=0;
                        pO_N = 0;
                        xO_N = 0;
                        etaO_N = 0;
                    else
                        pO_N = thetaO_N/((1+thetaO_N^ell)^(1/ell));
                        qO_N = pO_N/thetaO_N;
                        xO_N = VN(ind_phi,ind_h) - (c_effect/qO_N);
                        etaO_N = (xO_N - U(ind_phi,ind_h))/(VN(ind_phi,ind_h) - U(ind_phi,ind_h));
                    end  
                    t_VOhat_N = lambda*pO_N*xO_N + (1-lambda*pO_N)*VO(ind_phi,ind_h);
        
                    % For those Currently in N - Now Consider which job type is highest
                    if (t_VNhat_N>t_VNhat_O)
                        GN(ind_phi,ind_h) = 1;
                        opt_p_N(ind_phi,ind_h) = pN_N;
                        VN_hat_update(ind_phi,ind_h) = t_VNhat_N;
                        opt_x_N(ind_phi,ind_h) = xN_N;
                        opt_eta_N(ind_phi,ind_h) = etaN_N;
                        opt_theta_N(ind_phi,ind_h) = thetaN_N;
                    elseif (t_VNhat_O>t_VNhat_N)
                        GN(ind_phi,ind_h) = 2;
                        opt_p_N(ind_phi,ind_h) = 1;
                        VN_hat_update(ind_phi,ind_h) = t_VNhat_O;
                        opt_x_N(ind_phi,ind_h) = xO;
                        opt_eta_N(ind_phi,ind_h) = etaO;
                        opt_theta_N(ind_phi,ind_h) = 1.0;
                    else
                        GN(ind_phi,ind_h) = 3;
                        opt_p_N(ind_phi,ind_h) = pN_N;
                        VN_hat_update(ind_phi,ind_h) = t_VNhat_N;
                        opt_x_N(ind_phi,ind_h) = xN_N;
                        opt_eta_N(ind_phi,ind_h) = etaN_N;
                        opt_theta_N(ind_phi,ind_h) = thetaN_N;
                    end
        
                    % For those Currently in O - Now Consider which job type is highest
                    if (t_VOhat_N>t_VOhat_O)
                        GO(ind_phi,ind_h) = 1;
                        opt_p_O(ind_phi,ind_h) = pO_N;
                        VO_hat_update(ind_phi,ind_h) = t_VOhat_N;
                        opt_x_O(ind_phi,ind_h) = xO_N;
                        opt_eta_O(ind_phi,ind_h) = etaO_N;
                        opt_theta_O(ind_phi,ind_h) = thetaO_N;
                    elseif (t_VOhat_O>t_VOhat_N)
                        GO(ind_phi,ind_h) = 2;
                        opt_p_O(ind_phi,ind_h) = 1;
                        VO_hat_update(ind_phi,ind_h) = t_VOhat_O;
                        opt_x_O(ind_phi,ind_h) = xO;
                        opt_eta_O(ind_phi,ind_h) = etaO;
                        opt_theta_O(ind_phi,ind_h) = 1.0;
                    else 
                        GO(ind_phi,ind_h) = 3;
                        opt_p_O(ind_phi,ind_h) = pO_N;
                        VO_hat_update(ind_phi,ind_h) = t_VOhat_N;
                        opt_x_O(ind_phi,ind_h) = xO_N;
                        opt_eta_O(ind_phi,ind_h) = etaO_N;
                        opt_theta_O(ind_phi,ind_h) = thetaO_N;
                    end
        
                end % End loop over h values
            end % End loop over phi values
        
            error_U_hat = max(max(abs(U_hat_update(:,:) - U_hat(:,:))))/max(max(abs(U_hat_update(:,:))));
            error_VN_hat = max(max(abs(VN_hat_update(:,:) - VN_hat(:,:))))/max(max(abs(VN_hat_update(:,:))));
            error_VO_hat = max(max(abs(VO_hat_update(:,:) - VO_hat(:,:))))/max(max(abs(VO_hat_update(:,:))));
        
            vfi_error = max([error_U_hat,error_VN_hat,error_VO_hat]);
        
            
            U_hat(:,:) = U_hat_update(:,:);
            VN_hat(:,:) = VN_hat_update(:,:);
            VO_hat(:,:) = VO_hat_update(:,:);
            
        end % End while (vfi_error > vfi_toler)  
    
    
        % ================ Now that we have correct value functions, look at the expected value of entering O or N job types %

        for ind_phi = 1:nphi          % Loop over phi types
            for ind_h = 1:nh          % Loop over current h values
                c_effect = c;                                % Search cost
                    
                %---------------------------------------------------%
                % ------- Solve Search Decision of Unemployed ----- %
                xO = VO(ind_phi,ind_h) - c_effect;
                t_Uhat_O_compare(ind_phi,ind_h,ind_c) = xO;
        
                thetaN = (((VN(ind_phi,ind_h) - U(ind_phi,ind_h))/c_effect)^(ell/(1+ell)) - 1)^(1/ell);
                if (~isreal(thetaN))
                    pN = 0;
                    xN = 0;
                else
                    pN = thetaN/((1+thetaN^ell)^(1/ell));
                    qN = pN/thetaN;
                    xN = VN(ind_phi,ind_h) - (c_effect/qN);
                end
                t_Uhat_N_compare(ind_phi,ind_h,ind_c) = pN*xN + (1-pN)*U(ind_phi,ind_h);
       
            end % End loop over h values
        end % End loop over phi values


        
    end % end loop over c types 

end % End function